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component of particle velocity at the interface between 
media having different sound speeds and densities. 
Interface conditions are preserved for horizontal and 
sloping interfaces along a user-specified bottom profile. 


Test cases are included to demonstrate the use of the model. 
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LI. INTRODUCTION 


Since its introduction to the underwater acoustics 
community (Hardin and Tappert, 1973), the parabolic wave 
equation has stimulated a considerable amount of interest. 
The first solution programs used a split-step fast Fourier 
transform method to solve the parabolic equation; however, 
other solution techniques have been developed (McDaniel and 
Lee, 1982). One of the motives for developing alternative 
solution techniques is that problems arise when the Fourier 
transform encounters an interface between two media having 
different sound speeds or densities (Lee and Botseas, 1982). 

One alternative solution technique uses an implicit 
finite-~difference (IFD) solution method. The IFD method is 
unconditionally stable and has the capability to incorporate 
desired interface conditions. Implicit finite~difference 
methods for solving parabolic equations have been studied 
extensively by many authors. 

A computer program that utilizes the IFD method to solve 
the parabolic equation has been developed and is examined in 
detail in this thesis. The computer program predicts 
acoustic propagation loss in environments having user- 
specified bottom profiles. The program preserves continuity 


of pressure and continuity of the normal component of 





particle velocity at an interface between media having 
different sound speeds and densities. 

The program utilizes concepts developed by earlier 
authors. The use of the IFD method to solve the parabolic 
equation in underwater acoustics was developed by Lee and 
Papadakis (1979). The mathematical treatment of horizontal 
and sloping interfaces was developed by McDaniel and Lee 
(1982) and Lee and McDaniel (1983). And finally, the 
program utilizes some design features of an earlier computer 


program developed by Lee and Botseas (1982). 
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Il. PARABOLIC EQUATION 


A. INTRODUCTION 

The parabolic equation is an approximation to the 
elliptical wave equation. The derivation of the parabolic 
equation begins with the reduced wave equation (Helmholtz 


equation) in the form 


y*p+k*=0 De 
or 
ae p + acne = 1 Ze 
where 
k = wave number (= w/c) 
k, = reference wave number (= w/c,) 
n = index of refraction (= eye) 
p = time independent factor of complex pressure 
c = sound speed 
Co = reference sound speed 
w = angular source frequency (= 27 f£) 


For the case of cylindrical symmetry (2.2) becomes 
Cee Grae. © Dor + ka-n“p =U Des 
It is then assumed that p is of the form 
Paulo 2s (5) 
where u is a function of both range and depth and S is a 
DMP etoVmommnrangesoniy., Substitution of (2.4) into (2.3) 


and separation of variables shows that S(r) must satisfy 
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Bessel's equation of zero-order. For exp (-iwt) time 
dependence and outgoing waves, the solution is the zeroth- 
order Hankel function of the first kind, 

SiG) = Ho!) (kor). 
Further, u(r,z) must satisfy 


1 2 
+ (- + - § 


+ko(n¢ - 1) us0. 2.5 
E S 


u + U 


rr u 


ZZ x? ha 


With the help of the far-field asymptotic approximation for 
the Hankel function, (2.5) can be reduced to 


Gee, se 2ikou + ko“(n“ = 1) u = 0. 2.6 


Ei ZZ 


We now assume that u varies slowly with respect to range, 














ue | << +] 2 ae - eT 
combining (2.6) and (2.7) results in 
eeeeteeue = ko -(n- = 1) u = 0. 2.8 


Rearranging (2.8) results in the parabolic equation in the 


form 


uy = a(ky,r,z) u + b(k,l,F,z) u Phe, 


b ZZ 


where 


deeryiz) = (iko/2) Imo(e,c)o= 11 


b(k,,r,z) = i/2k,. 


The assumption (2.7), fundamental to the parabolic equation 


method, is equivalent to neglecting back-scatter. 


B. SPLIT-STEP FAST FOURIER TRANSFORM SOLUTION 


1. Description 


For the first few years after its introduction into 


the acoustical community the parabolic equation was solved 


WZ 





exclusively with the help of the split-step fast Fourier 
transform (SSFFT) method developed by Tappert and Hardin 
(gensen and Krol, 1975). In this method, UD eee 8 ) ais 
represented by the inverse transform of its Fourier 
transform. The SSFFT method requires periodic boundary 
conditions in z because of the finite Fourier transform. 
iiss is handled bYyointroducing an artificial, horizontal, 
pressure release bottom below the physical bottom. The 
SoPrT method is umconditionally stable (Brock, 1978). 

The SSFFT method has been implemented by Jensen and 
Krol and by Brock. Detailed descriptions can be found in 
publications of Jensen and Krol (1975) and Brock (1978). 

2. Interface Treatment 

paces fiemte roaucea by “the SSFFT method are 
proportional to the range step and to eee is Eire 
index of refraction (Jensen and Krol, 1975). Because the 
index of refraction has a large change across the water- 


sediment interface n will be large and thus the error will 


Ze, 
bee anges To reduce this error, a very small horizontal 
range step must be used. However, this results in very long 
computer execution time. The problem of a discontinuity in 
sound speed at the water-sediment interface and the 
resultant difficulties in solving shallow water propagation 


problems using the SSFFT method are addressed in Jensen and 


Krol (1975). 
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Another more serious problem with the SSFFT method 
is that it neglects any density difference between the water 
and the sediment. A density difference can be important in 
that it influences the reflection coefficient. The problem 
becomes more significant as the density difference becomes 
larger. 

In summary, the discontinuities in sound speed and 
in density at the water-sediment interface cause problems 
for the SSFFT method. The SSFFT method is therefore 
intrinsically better suited for deep water propagation 
environments for which the water-sediment interface is an 


unimportant feature. 


C. IMPLICIT FINITE-DIFFERENCE SOLUTION METHOD 

In 1979 Lee and Papadakis introduced the Crank-Nicolson 
implicit finite-difference method to solve the parabolic 
equation for underwater acoustic propagation. The Crank- 
Nicolson method uses a second-order central difference 
formula to approximate u,, in (2.9) and casts the problem in 
the form of a tridiagonal matrix system. A representative 


th 


row in the matrix system (the m row) is 
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m=-1 
eke ntl Be rie AGG 2). TT) ae ot 7 | n+l 
( 2 n2 Onn » 1 - 2 Hs orn i 2 a ee n2 DT ) m 
ntl Pee G 
m+1 | 
un 
m—-1 
lode n i kK n 1k n n 
==-—., b, eae cede =. =-—-, bb) u 
Zz h2 m 2 h2 m 2 h2 m i 
Yat 


where 

k = horizontal range increment 

h = vertical depth increment 
Supercripts are used to indicate range indices and 
Subscripts are used to indicate depth indices. In (2.10) 
the field is known at range index n and is to be solved at 
range index n+ 1. Therefore, the right hand side of (2.10) 
reduces to a single, known value and the solution field is 
advanced from range index n to range index n +1 by solving 
the tridiagonal system of equations. 

The IFD scheme is consistant, unconditionally stable and 
it converges to the theoretical solution as the range and 
depth increments tend to zero (Lee et al., 1981). An 
advantage of selecting an implicit scheme over an explicit 
scheme is that an explicit scheme is only conditionally 
stable (Lee and Papadakis, 1979). Another advantage of the 


implicit scheme is smaller errors. More detailed 


iss 





information addressing both implicit and explicit solutions 
of parabolic equations can be found in Gerald (1980). 

The first IFD scheme handled discontinuities in the 
speed of sound profile but did not consider the effects of 
density discontinuities. It therefore did not correctly 
treat the interface between media having different 
densities. 

D. IMPLICIT FINITE-DIFFERENCE METHOD: TREATMENT OF A 

HORIZONTAL INTERFACE 

In 1982 McDaniel and Lee introduced a scheme for 
micOmpOoratbing:a eeticoa tal interface into the IFD method. 


The interface separates two media with different sound 


speeds and densities (Figure 1). 


Z tis 
AN 


Figure 1. IFD Treatment of a Horizontal Interface 
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The scheme preserves continuity of pressure and continuity 
of the normal component of particle velocity across the 
interface and does not affect the stability of the IFD 
method. 

The interface equation in the tridiagonal matrix system 
that results from incorporating continuity of pressure and 
continuity of the normal component of particle velocity is 


(McDaniel and Lee, 1982) 


yar 
(x pntl- 1-2, pati gntl a 
h2 m Z m m 
yitl 
m 
‘ mo “Sh, - - ped ntl 
h LA, h 2 He 


where 4 
Py 1 
po | a 
by F252 
ay 1 42 
OS OY i el 
b, -2 b2 
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ae 
a, 


a and b are defined in (2.9) 


density in layer 1 (water) 


density in layer 2 (sediment) 


Incorporating the horizontal interface into the IFD method 
requires inserting (2.10) for the row in the tridiagonal 
matrix system that corresponds to the interface. 

The error in the solution is 

oO (k? + kh) 
on the interface and 

O (k? + kh?) 
in a continuous medium (McDaniel and Lee, 1982). 

E. IMPLICIT FINITE-DIFFERENCE METHOD: TREATMENT OF A 

SLOPING INTERFACE 

In 1983, Lee and McDaniel extended their treatment of an 
interface between two media to include the case of a sloping 
interface. As for the case of a horizontal interface, the 
treatment of a sloping interface preserves continuity of 
pressure and continuity of the normal component of particle 
velocity at the interface between media having different 
densities and sound speeds. 

The problem of a sloping interface is separated into two 
cases: downslope, and upslope. Each case requires inserting 
two new rows into the IFD tridiagonal matrix system. One 
new row is required at the level corresponding to the 
interface at the range where the solution is known and the 


second new row is required at the level corresponding to the 
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interface at the range at which the solution is to be 
solved. Therefore, four new equations are required to 
cover both the downslope and upslope cases. The four 
sloping interface equations are derived and shown in Lee and 
McDaniel (1983). The equations are somewhat involved but 
they are of the same tridiagonal form as the original IFD 
matrix equations. The error for points on or adjacent to a 
sloping interface is (Lee and McDaniel, 1983) 


O (k? + kh). 


1.9 





III. COMPUTER IMPLEMENTATION 


A. INTRODUCTION 

An IFD solution program that implements (2.10), (2.11) 
and the four sloping interface equations has been developed 
to solve underwater propagation problems. The program is 
written in FORTRAN using single precision, complex 
arithmetic and has been installed on the Naval Postgraduate 
School's IBM-3033 digital computer. Appendix A contains a 
program listing. 

The solution program consists of one main program and 20 
subroutines. A modular program construction was selected 
Eemetlexibi1lity and clarity. 

Within each routine a structured programming approach is 
Meal 12 ed. The structured program format, coupled with 
generous commenting, makes the program relatively easy to 
trace through. 

As presently installed on the NPS computer the solution 


program is run interactively. Appendix B contains details. 


B. GENERAL CHARACTERISTICS 

The solution program handles the following environmental 
conditions: range independent sound speed profile in the 
water column, range dependent bottom profile and iso-speed, 


iso-density sedimentary bottom layer. The program utilizes 
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a Gaussian starting field and an artificial pressure release 
surface in the sediment at a user-specified depth. (An 
artificial pressure release surface is not required for the 
IFD solution method; however, such a surface permits 
efficient solution for the pressure field when it is known 
to tend to zero at great depths in the bottom.) 

Attenuation in the water and sediment is introduced 
using complex indices of refraction. Artificially strong 
attenuation is applied in the lower portion of the sediment 
layer to enhance attenuation of the field above the 
artificial pressure release surface. 

A user-specified bottom profile is input as a series of 
linear bottom segments. The range step along a horizontal 
bottom segment is set to a user-specified value. The range 
step along a sloping bottom segment is automatically set by 
the program so that the sloping bottom intersects the next 
vertical grid point. Bottom modifications are required in 
certain situations to meet the requirements that the range 
step not be too large and that the interface must pass 
through a grid point at every range at which the pressure 
field is solved. For avery gently sloping bottom, if the 
calculated range step exceeds a user-specified value then 
the program will automatically model the bottom as a series 
of level and sloping sections. For these cases, the 
difference between the modified bottom and the user- 


specified bottom is always less than or equal to one-half 
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the vertical grid increment. The user is informed if the 
bottom is modified. 

It is foreseen that future enhancements will increase 
the program's generality. In particular, changes to allow 
range dependent sound speed profiles, sound speed profiles 
in the sediment and a user-specified starting field should 
be relatively simple. The modular construction of the 


program facilitates these types of changes. 


C. MAIN PROGRAM IFD 

IFD is the main program. It controls program execution 
and calls subroutines as appropriate. 

The first executable statement in IFD calls subroutine 
ERRSET, a system subroutine peculiar to the NPS computer 
that correctly sets a variable value to zero when an 
underflow condition exists. Most computer systems do this 
automatically; however, depending on the particular system a 
call similar to ERRSET may be required. 

IFD then calls subroutine READ to read input data, 
subroutine SVPW to calculate the sound speed at grid points 
in the water column, subroutine INITAL to initialize 
constants and variables, and then subroutine MATCON to 
calculate matrix constants. Subroutine SFIELD is then 
called to calculate the Gaussian starting field, followed by 
subroutines WRITE1 and PRINT1 which write and print output 


data. In the context of this program, "write" refers to 


ee 





writing unformatted data into a file to be used by the 
plotting routine and "print" refers to writing formatted 
data into a file which can be sent directly to the printer. 
IFD then calls subroutine NEWSEG which is the beginning 
of a loop that is called every time a new linear bottom 
segment is reached. NEWSEG calculates variables that 
characterize a new bottom segment. The next call is to 
subroutine NEWMAT which calculates matrix elements for the 
new bottom segment and advances the solution field one range 
step. IFD then enters a loop that advances the solution one 
range step for every pass through the loop. Inside the loop 
the range markers are advanced and the solution is advanced 
one step for the downslope, level, upslope, modified bottom 
downslope or modified bottom upslope situation as 
appropriate. In addition, the artificial attenuation 
mentioned earlier is applied by calling ATTENU and calls are 
made to WRITE2 or PRINT2 as required. Finally the range is 
checked to see if it has advanced to the maximum range 
specified for the problem. If it has, then IFD calls 
subroutine END which passes appropriate messages to the 


terminal and stops program execution. 


D. SUBROUTINES 
1. Subroutine READ 
Subroutine READ is called by IFD to read input data 


from unit number NIU = 51. Input data are read in free 
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format and data are transferred back to main program IFD via 
common blocks. READ contains error checks for (1) input 
data insufficient and (2) the final depth in the sound speed 
profile unequal to the maximum depth in the water column. 
If either of these error conditions exists, READ issues an 
appropriate error message to the terminal and stops 
execution, 
2. Subroutine SVPW 
Subroutine SVPW calculates the vertical grid spacing 
used in the water and sediment. It also calculates the 
speed of sound at each of the grid points in the water 
column. Linear interpolation is used to calculate the sound 
speed at grid points between points on the user-specified 
sound speed profile. 
3. Subroutine INITAL 
Subroutine INITAL initializes constants and 
variables. If the user inputs 0.0 for the value of the 
reference sound speed then INITAL sets the reference sound 


speed c. to the sound speed averaged over the deepest water 


O 
column. If the user inputs 0.0 for the value of the maximum 
range step then INITAL sets the maximum range step to the 
reference wavelength, 

DRMAX = XLAMDA . 
Setting the maximum range step to the reference wavelength 


is somewhat arbitrary; however, until the actual limit on 


the range step is better understood it serves as a rough 
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mere OL (ehumb, Finally, if the user inputs 0.0 for the 
value of the range step along a horizontal interface then 
INITAL sets the range step to half of the reference 
wavelength, 
DRLVL = 0.5 * XLAMDA. 
The default range step along a horizontal interface is half 
the default maximum range step. 
4, Subroutine MATCON 

Subroutine MATCON calculates constants needed to 
compute tridiagonal matrix elements. Most of the constants 
computed in MATCON have no direct physical significance but 
contribute to computational efficiency. Attenuation in both 
the water and sediment is calculated with the help of a 


complex index of refraction n, 


on BETA 
22) |) ey eae 
a 54575054 
or 
2 Z 
> So . Co BETA 
n=: jroo + 1 — es 
C+ Cy Zee OD 
where 
BETA = attenuation (dB/wavelength) 
ce = reference sound speed (m/s) 
Cy = sound speed (m/s) at point j 
54.575054 = conversion factor used in converting 


db/wavelength to nepers/meter. 
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5. Subroutine SFIELD 
Subroutine SFIELD calculates the Gaussian starting 
field at range r = 0. This subroutine is identical with 
that of. Lee and Botseas (1982); both yield the starting 


field suggested by Brock (1978), 


U(I) = CMPLX (PR, 0.0) 
where 
ZM-zs | 2 -~ZM-ZS |? 
GW GW 
PR = GA [e -e ] 


ZM = depth (m) of grid point 

ZS = source depth (m) 

GW = Gaussian width (m) (= 2/FK) 

FK = reference wave number (1/m) (= 2% f/c,) 


GA = Gaussian amplitude [= (2/FK)!/2/cew 


6. Subroutine WRITE! 

Subroutine WRITE1 outputs data toa file that is used 
by the plotting routine. This output file corresponds to 
unit file number NOU = 52. 

eoGoucIne PRINT | 

Subroutine PRINT1 outputs formatted data toa file 

that can be sent to the printer. Ene Output ft le 


corresponds to unit file number NPOUT = 55. 
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8. Subroutine NEWSEG 

Subroutine NEWSEG is called at the start of each new 
linear bottom segment. NEWSEG computes and initializes 
variables that depend on characteristics of the segment. 
One of the variables initialized is ISLOPE which is a slope 
flag having value 1 if the bottom slopes down, 2 if the 
beoeeom 1s horizontal, 3 if the bottom slopes up, 4 if the 
bottom slopes down but must be modified because the slope is 
Vemy small, or 5S if the bottom slopes up but must be 
modified because the slope is very small. NEWSEG also 
issues error or warning messages as appropriate. 

9. Subroutine NEWMAT 

Subroutine NEWMAT calculates matrix elements for the 
X and Y matrices. The Y matrix corresponds to the range at 
which the solution field is known and the X matrix 
corresponds to the range at which the solution field is to 
be found. The Y matrix is multiplied by the known solution 
field to obtain the right-hand side column vector needed to 
solve the tridiagonal system. 

NEWMAT sets up the tridiagonal matrix system for the 
new bottom segment and than calls TRIDG to solve the system 
at the first range step. It then calls ATTENU to apply 
artificial attenuation, calls WRITE2 or PRINT2 as required, 
and finally updates the interface pointer that indicates the 


index of the grid point at the water-sediment interface. 
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10. Subroutine WRITE2 
Subroutine WRITE2 is basically a continuation of 
subroutine WRITE]. It outputs data to a file corresponding 
to unit file number NOU = 52. This is the file that is 
used by the plotting routine. At range intervals specified 
by the user, WRITE2 outputs range, depth, and the value of 
mer, Z ) « 
11. Subroutine PRINT2 
Subroutine PRINT2 is basically a continuation of 
subroutine PRINT1. It outputs formatted data to a file 
corresponding to unit file number NPOUT = 55. This file can 
be sent directly to the printer. At range and depth 
intervals specified by the user PRINT2 outputs tabular 
values of transmission loss and u(r,z). 
12. Subroutine TRIDG 
Subroutine TRIDG solves a linear tridiagonal matrix 
system. TRIDG is a modified version of subroutine TRID as 
listed in Gerald (1980). The major modifications to 
subroutine TRID involved introducing arrays CTWO and CR to 
preserve the original matrix element values and to make the 
routine more efficient. Introducing the two new arrays 
requires more storage space but results in a substantial 
savings in execution time, particularly for the case of a 


horizontal interface. 
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13. Subroutine TRIDL 
Subroutine TRIDL is a modified version of subroutine 
TRIDG,. TRIDL differs from TRIDG in that it does not compute 
arrays CTWO and CR but rather uses the array values 
calculated in TRIDG. 
14. Subroutine DOWN 
Subroutine DOWN updates the tridiagonal matrix and 
calls subroutine RHS to update the right-hand side for the 
case of a downward sloping interface. DOWN then calls 
subroutine TRIDG to solve the tridiagonal system of 
equations and finally updates the interface pointer. 
15. Subroutine UP 
Subroutine UP performs exactly the same tasks as 
subroutine DOWN, but for the case of an upward sloping 
interface. Subroutine TRIDG is again called to solve the 
system. 
16. Subroutine LEVEL 
Subroutine LEVEL is called to advance the solution 
for the case of a horizontal interface. For this case the 
tridiagonal matrix elements at the advanced range need not 
be changed from the previous calculation. Therefore, LEVEL 
need only update the right-hand side by calling RHS and 
then solve the system by calling TRIDL. 
17. Subroutine RHS 
Subroutine RHS computes the right-hand side of the 


tridiagonal system by multiplying tridiagonal matrix Y by 
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the known solution field U(I). The resultant right-hand 
side column vector is stored in C(I,4). 
18. Subroutine SSLOPE 
Subroutine SSLOPE is called to advance the solution 
in the case where the bottom has been modified. SSLOPE 
determines which of three cases a particular section falls 
into: a level section following a level section, a level 
section following a sloping section, or a sloping section. 
For the case of a level section following a level section 
SSLOPE calls LEVEL to advance the solution. For the case of 
a level section following a sloping section SSLOPE updates 
appropriate matrix elements, calls RHS and then calls TRIDG 
to advance the solution. And in the case of a sloping 
section SSLOPE updates matrix elements and calls either DOWN 
or UP as appropriate. 
19. Subroutine ATTENU 
Subroutine ATTENU applies artificial attenuation to 
the bottom portion of the sediment layer as suggested by 
Brock (1978). The artificial attenuation matrix ATT(1) is 
calculated in subroutine NEWMAT. 
20. Subroutine END 
Subroutine END is called when the solution field has 
reached the maximum range specified. END sends appropriate 
messages to the terminal and stops execution. The messages 


are applicable to the program as installed on the NPS 
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computer but may not be appropriate for the program if 


installed on another system. 


E. INPUT DATA 
1, Input File 
The input data must be stored in a file 
corresponding to unit number NIU as assigned in subroutine 
READ. In its present form READ sets the input unit number 
to NIU = 51. If the user prefers to read the data from a 
different unit (for example, a card reader), then variable 
NIU in READ should be set equal to the appropriate unit 
number. 
Ze lnpae Format 
The input data is read in free format. The input 
card images (or input cards) are arranged as follows: 


CARD CONTENTS 


1 EROs.2S, ZR, CO, N 
where 
FRQ = frequency (Hz) 
ZS = source depth (m) 
ZR = receiver depth (m) 


(program will reset to depth of nearest 
grid point) 

CO = reference sound speed (m/s) 
If cO = 0.0, CO is set to the sound speed 


averaged over the deepest water column. 
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number of vertical grid points 


If the user desires that every integer 


depth value correspond to a grid point then 


(neglecting dimensions) N should be set to 


an integer multiple of ZLYR2, the depth of 


the pressure release surface. 


CARD CONTENTS 


2 RMAX, DRLVL, DRMAX, WDR, PDR, PDZ 
where 

RMAX = maximum range (m) of solution 

DRLVL = range step (m) for marching solution 
along horizontal interface 
If DRLVL = 0.0, then DRLVL is set to 1/2 
wavelength. 
If DRLVL is greater than DRMAX, then 
DRLVL is set to DRMAX, 

DRMAX = maximum allowable range step (m) 
Tf DRMAX = 0.0, then DRMAX is set to 1 
wavelength. 

WDR = range increment (m) at which solution is 
written to file used by plotting routine 

PDR = range increment (m) at which solution is 
printed 

PDZ = depth increment (m) rounded to nearest DZ 


at which solution is printed 
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CARD 


N+1 


N+2 


N+3 


N+4 


N+M 


CONTENTS 


Bayt), B2Z(1) 


BR(2) , BZ(2) BOTTOM PROFILE 


BR(3) , BZ(3) Range and depth of water (m). 


Maximum number of points = 100, 
Program will reset depths to 


nearest grid point. 


-—1 This card marks end of bottom 


profile. 


ZLYR1, RHO1, BETA1 


where 


ZLYR1 


RHO1 


BETA1 


ZSVP(1), 


ZSVP(2), 


maximum water depth (m) 

density of water (g/cm?) 

attenuation in water (dB/meter) 

If BETA1 is less than 0.0 then program 
calculates BETA1 with an empirical 
formula (Brock, 1978). 


csvP(1)] SOUND SPEED PROFILE 


CSVP(2) Depth (m) and sound speed 
‘ (m/s). 
: ZSVP(1) must equal 0. 


The last depth must equal 


ZLYR1. 
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CARD CONTENTS 


N+M+1  ZLYR2, RHO2, BETA2, C2 


where 
ZLYR2 = depth (m) of pressure release surface at 
bottom of sediment layer 
RHO2 = density (g/cm?) of sediment 
BETA2Z = attenuation (dB/wavelength) in sediment 
eZ = sound speed (m/s) in sediment 


N+M+2 ZABLYR 
where 


ZABLYR 


depth (m) of upper surface of artificial 
attenuation layer in sediment. 


ZABLYR should be about 3/4 of ZLYR2. 


F. PROGRAM OUTPUT 
1. Output Printer File 
The program outputs formatted data to a file 
corresponding to unit number NPOUT which is set to 55. This 
formatted data file may be sent to the printer if desired. 
On another system the user may elect to assign NPOUT to the 
unit number corresponding to the printer and thereby send 
the formatted data directly to the printer. 
Zoeeourout Pletter File 
The program outputs unformatted data to a file that 


is used by the plotting routine. The unit number for this 
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file is NOU which is set to 52. The output data in this 


file are stored as follows: 


LINE CONTENTS 


4 RMAX 
where 
RMAX = Maximum range (m) of solution 
Z rs, Lik eu 
3 RA, ZR, U 
where 
RA = range (m) 
ZR = depth (m) 
U = complex u at specified range and depth 


3. Terminal Output 
Certain WRITE statements in the program specify unit 
Number 6. Unit number 6 on the NPS computer for an 
interactive program corresponds to terminal output. If the 
program is not run interactively then WRITE statements with 
unit number 6 may be deleted. Any pertinent information 
passed to the terminal during an interactive run is also 


passed to unit number 55. 


G. PLOTTING PROGRAM 
Appendix C contains a listing of the IFD plotting 


program installed on the NPS computer. The filename and 
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filetype of the program are PLOTIFD FORTRAN. This program 
was written separately from the IFD program because it was 
recognized that different computer installations have 
different plotting facilities. For users with different 
facilities the program will be a helpful reference. 

The program reads data from unit number NOU = 52 which 
corresponds to a file with filename/filetype IFDOUT PLOTTER. 
Details concerning using the plotting routine are included 


in Appendix B. 
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LV. PEST CASES 


A. HORIZONTAL INTERFACE CASES 

The IFD program treats a horizontal interface using the 
same theoretical approach as the IFD program published by 
Lee and Botseas (1982). Throughout the remainder of this 
thesis the Lee and Botseas (1982) program will be called the 
FINITE-DIFFERENCE program and the program presented in this 
thesis will be called the IFD program. Two cases were run 
to confirm that the IFD program is in agreement with the 
FINITE-DIFFERENCE program for horizontal interfaces. 

1. Isospeed Shallow Water 

This case, first published by Jensen and Kuperman 

(1979), considers propagation in a shallow water, isospeed 
environment. The water depth is 100 meters and the solution 
field is calculated out to 25 kilometers. The sound speed 
is 1500 m/s in the water and 1550 m/s in the sediment. 


3 and 1 


Density and attenuation in the sediment are 1.2 g/cm 
dB/wavelength respectively. The source and receiver are 
both at 50 m and the source frequency is 500 Hz. 

Solutions obtained uSing a normal mode program 
(SNAP), a split-step fast Fourier transform program (PAREQ) 
and the FINITE-DIFFERENCE program are shown in Figure 2. 
SNAP and PAREQ are programs that were developed at SACLANT 


Centre and are discussed further in Jensen and Kuperman 
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Figitwre 2. 
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no? 9). The solution obtained using the IFD program is 
shown in Figure 3. All of the solutions are in excellent 
agreement. 

The input runstream that produced the results shown 
in Figure 3 for the IFD program is as follows: 


Input Runstream 


500 50 50 0 500 

25000 5 » 50 5000 50 
0 100 

25000 100 

-1 -1 

100 dd =120 

0 1500 

100 1500 

20 VS uae) 150 

200 


Zaeorizontal” Interface 

This case, called the "horizontal interface problem" 
in Lee and Botseas (1982), considers propagation in an 
environment with the sound speed profile shown in Figure 4. 
Source frequency is 100 Hz, source depth is 30 m and receiv- 
er depth is 90 m. The density in the bottom is 2.1 g/cm? 
and the sound speed in the bottom is 1505 m/s. No attenua- 
tion is applied in the water or sediment using complex 


indices of refraction; however, artificial attenuation is 


applied in the lower portion of the sediment layer. 
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SOUND SPEED (m/s) 
1480 1490 1500 1510 
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FREQUENCY = 100 Hz 


SOURCE DEPTH = 30m 
480 


512 RECEIVER DEPTH = 90m 


Figure 4. Horizontal Interface Case 


The solution obtained using the IFD program is shown 
in Figure 5. This solution is in excellent agreement with 


the FINITE-DIFFERENCE solution shown in Lee and Botseas 


(loZ).. 
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The input runstream for this case is as follows: 


Input Runstream 


100 30 90 0 600 
20000 Z Z 50 10000 50 
0 240 

20000 240 

-1 -1 

240 1 18, 0.0 

0 1500 

A0, 1498 

240 1500 

1200 Za 0.0 150 
Se 


Bee RANGE-DEPENDENT CASES 

The following range-dependent cases were solved by 
Jensen and Kuperman using SNAP, the normal mode program, and 
PAREQ, the SSFFT program (Jensen and Kuperman, 1979). The 
cases were also solved by Lee and Botseas using the FINITE- 
DIFFERENCE program (Lee and Botseas, 1982). The FINITE- 
DIFFERENCE program treats the sloping interface as a "Stair 
step" and uses the interface conditions appropriate for a 
horizontal interface. The IFD program handles the sloping 
interface using the interface treatment derived by Lee and 


Memaniei (1983). 
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1. Deep-to-Shallow Water 
This case considers propagation in an environment as 
@epeeted in Bigure 6. The problem is solved for a bottom 
with a 8.5 degree upslope and one with a 0.85 degree 
upslope. source frequency is 25 Hz, source depth and 
meee Zem depen are 25 meters. Sound speed in the water is 
1500 m/s. In the sediment, sound speed is 1600 m/s, density 


moe 1.5 g/em? and attenuation is 0.2 dB/wavelength. 


O lO 20 2© 40 km 






= 50 siiciesteacteiieietats 

= oop C2 1500 m/s 

= ie Cy 2 1600 m/s 

T 250 P5158 Yom 

CO 300 B+ 0.2 d8/A 
390 


Figure 5. Deep-to-Shallow Water Case 


The results for the 8.5 degree upslope case as 
produced by SNAP,PAREQ, and FINITE-DIFFERENCE are shown in 
Figure 7. The results as produced by IFD are shown in 
Figure 8. The difference between the results produced by 
SleeeancG PAREO is attributed to failure of the “adiabatic” 
theory underlying the SNAP program (Jensen and Kuperman, 
1979). As determined from the input runstream the FINITE- 
DIFFERENCE results were obtained using 1.0 rather than 1.5 
g/em? For the density in the sediment (Lee and Botseas, 


Po 2 4. 


a4 





5OQ 


60 


ie 


80 


90 


100 


110 


PROPAGATION LOSS (dB) 


120 


130 


SO 


60 


70 


80 


90 


100 


PROPAGATION LOSS (dB) 


L2e@ 


130 


~~ = oe * 





*s Gd = 8.5° 


10 20 30 40 
RANGE (km) 


FINITE—OIFFERENCE 


10 15 20 25 30 35 40 


RANGE (km) 


Figure 7. Propagation Loss Versus Range for Deep-to-Shallow 
Water Case, 8.5 Degree Slope; SNAP, PAREO and 
FINITE-DIFFERENCE Results 
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Pump ue rumstream Lor the ELFD program that 
produced the results shown in Figure 8 is as follows: 


Input Runstream 


25 25 25 0 1000 
40000 10 0 100 10000 25 
0 350 

10000 350 

12000 50 

40000 50 

= St 

350 veye =f 

0 1500 

350 1500 

1000 1.5 ee 1600 

750 


Appendix D contains the printed output produced by the IFD 
program using the above runstream. 

The results for the 0.85 degree upslope case as 
produced by SNAP and PAREQ are shown in Figure 9. The 
results produced by IFD are shown in Figure 10. 

2. Shallow-to-Deep Water 

This case considers propagation in the environment 
depicted in Figure 11. The environment is exactly the same 
as the deep-to-shallow water environment except that the 
shallow and deep portions have been reversed and thus the 


bottom slopes down rather than up. 
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Figure 9. Propagation Loss Versus Range for Deep-to-Shallow 
Water Case, 0.85 Degree Slope; SNAP and PAREQ 
Results 
The results for the 8.5 degrees downslope case as 
produced by SNAP and PAREQ are shown in Figure 12. As 
before, the difference between SNAP and PAREO is attributed 
to failure of the SNAP program. The results produced by IFD 
are shown in Figure 13. 
The results for the 0.85 degree downslope case are 
shown in Figures 14 and 15. 
3. Comments 
Differences between the results obtained using the 
SNAP and PAREQ programs for the range-dependent cases are 
discussed in Jensen and Kuperman (1979). The major 
Gimme rmences are attributed to the violation cf the 


adiabatic assumption in the SNAP program. 
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Poegure 14, Propagation Loss Versus Range for Shallow-to- 


Deep Water Case, 0.85 degree Slope; SNAP and 

PAREQ Results 
Figure 16 shows the results produced by IFD for the 
8.5 degree, deep-to-shallow water case if 1.0 rather than 
es g/cem? is used for the density of the sediment. The IFD 
results obtained using 1.0 g/cm? are in very close agreement 
with the PAREQ results shown in Figure 7. However, when the 
correct value of 1.5 g/cm? is used the slope of the 
propagation loss curve beyond 15 km is less steep (Figure 8) 
and the results do not agree as well with the results 
produced by PAREQ. The observation that the slope of the 
propagation loss curve becomes less steep when 1.5 g/cm? is 


used is qualitatively consistant because a higher density 


difference means a higher reflection coefficient which in 
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turn results in more energy confined within the water layer. 
The difference in slope between the PAREQ and IFD results is 
attributed to the PAREQ program's apparent failure to 
account for the affects of the density discontinuity at the 
water-sediment interface. Because the IFD program correctly 
accounts for density discontinuities the results produced by 
IFD are believed to be more accurate than those of PAREQ 


when interface interaction is important. 
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V. COMMENTS AND CONCLUSIONS 


(re TER method 1S an efficient, stable method for 
solving the parabolic equation. Use of the IFD method is 
particularly advantageous in shallow water environments 
where the water-sediment interface is an important 
parameter. 

The IFD program presented in this thesis incorporates 
G@entimuity of pressure and continuity of the normal 
component of particle velocity across horizontal and sloping 
interfaces. The program's capability to incorporate the 
exact interface conditions on a sloping interface, to 
automatically determine step-size, and to modify the bottom 
as required for the case of avery gently sloping bottom are 
important features. 

Projected program enhancements include wide angle 
propagation (Lee and Gilbert, 1982), range-dependent sound 
speed profiles in the water, range-dependent sound speed 
profiles in the sediment layer, and multiple sediment layers 
with horizontal or sloping interfaces. These enhancements 
are listed in their approximate order of importance. The 
program's modular construction and structured style will 


facilitate implementation of these enhancements. 
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APPENDIX A: IMPLICIT FINITE-DIFFERENCE PROGRAM LISTING 
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APPENDIX B 


RUNNING THE IMPLICIT FINITE-DIFFERENCE 
PROGRAM ON THE NPS COMPUTER 


eee LRODUCT ION 
This appendix describes one simple procedure for running 


the IFD program on the NPS computer. 


B. COPYING FILES ONTO USER'S DISK 
Four files are needed to run the IFD program; the 


filenames and filetypes are 


1B, FORTRAN 
1B, BAEC 
PEOrTED FORTRAN 
EEOTEED EXEC 


The files are available from a computer account maintained 
by the underwater acoustics curriculum. To link with this 
account and obtain copies of the files the user should 
proceed as follows: 

(1) Log on terminal 

(2) Enter: CP LINK 0160P 191 195 RR 

(3) When prompted for read password enter: UX 

wie enter: ACC 195 C 

eeemcer:s COPYIFD 
At this point the four files should reside on the user's A 


dvsk. 
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C. RUNNING THE IFD PROGRAM 

Before running the IFD program the user should define 
additional storage space, compile the program, and set up 
the input data file. To define additional storage space, 


(1) Enter: DEF STOR 1M 
2) Senter: [£ CMS 


These two commands need only be entered one time for each 
terminal session; the storage space will remain for the 
entire session. 

To compile the program, 

Enter: FORTGI IFD 
The program need be compiled only one time unless the 
program is changed, in which case the new version should be 
recompiled. 

The final step before running the program is setting up 
the input data file which has filename and filetype IFDIN 
DATAIN. The user must create or modify this file so that it 
contains input data as described in Section III.E.2 of this 
thesis. For more information concerning how to create or 
modify files see NPS Technical Note TN-VM-02 which is 
available in the computer consultant's office. 

If the above steps are accomplished, the user can then 
run the program with 

Enter. Lr p 
Shortly after entering this command the user will be 


prompted for a run identification. The run identification 
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is an arbitrary identification label that will appear on the 
output printer file. Enter run identification as desired. 
At the end of the run the user is informed of the two 
output data files generated by the program. If desired by 
the user the output printer file (IFDOUT PRINTER) may be 
sent to the printer. The output plotter file, IFDOUT 


PLOTTER, serves aS an input file for the plotting program. 


D. RUNNING THE PLOTTING PROGRAM 

A Tektronix-618 terminal is used to run the plotting 
program. The first step is to log onto the terminal in the 
normal manner and then define additional storage space by, 


fie eeiter: DEF STOR 1M 
C7) ) emeer: I CMS 


The plotting program has filename and filetype PLOTIFD 

FORTRAN. To compile the program, 
Po@ems TORT GE PLOT fr D 
Unless the program is changed it need only be compiled one 
time. To run the program 
ghecia 3 [210i >). 

The user will be prompted for axes and smoothing 
information; enter values and responses as appropriate. The 
transmission loss curve will be displayed on the CRT screen, 
a hard copy may be obtained by pushing the HARD COPY button 
under the screen, and the screen may be cleared by pushing 


the ENTER key. 
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